REMARKS 



Applicant intends this response to be a complete response to the Examiner's 7 May 2009 
Non-Final Office Action. Applicant has labeled the paragraphs in his response to correspond to the 
paragraph labeling in the Office Action for the convenience of the Examiner. 

DETAILED ACTION 
Continued Examination Under 37 CFR 1.114 
The Examiner states as follows: 

1 . A request for continued examination under 37 CFR 1.114, including the fee set forth in 37 

CFR 1.1 7( e), was filed in this application after final rejection. Since this application is eligible for 
continued examination under 37 CFR 1 .114, and the fee set forth in 37 CFR 1 .1 7(e) has been timely 
paid, the finality of the previous Office action has been withdrawn pursuant to 37 CFR 1.114. 
Applicant's submission filed on March 23,2009 has been entered. 

Applicants acknowledge the Examiner's statements. 



Response to Amendment and Arguments 

The Examiner states as follows: 

2. Applicant's amendment filed on March 2.^1.2009 has been entered and made of record. 

Applicants acknowledge the Examiner's statements. 

Requirement for Information under 37 CFR 1.105 

The Examiner states as follows: 

3. Applicant and the assignee of this application are required under 37 CFR 1.105 to provide 
the following information that the examiner has determined is reasonablynecessary to the examination 

of this application. 

In response to this requirement, please provide a copy of the article "Nonseparable 
bidimensional wavelet bases" by Cohen et. al.. which was cited in the article "Non-separable Radial 
Frame Multiresolution Analysis in Multidimensions and Isotropic Fast Wavelet Algorithms" by 
Papadakis et al. In addition, please provide a concise explanation of the reliance placed on that 
publication in the development ofthe disclosed subject matter. 

Applicants points out to the Examiner that most of the cited references are not relevant to the 
invention, but merely direct the reader to where certain mathematic derivations and/or functions are 

described or where certain mathematical understandings can be found and where certain 
understandings about human vision can be found. The following list of references and there specific 
usage in the specification are set forth below: 



[6] J.J. Benedetto and S. Li. The Theory of Multiresolution Analysis Frames and Applications 
to FilterBanks. Appl. Comp. Harm. Anal., 5:389-427, 199&, deal with filters that can be constructed 
from the filters of this invention. 



[7] A. Cohen and I. Daubechies. Nonseparable Bidimensional Wavelet Bases. Revista 
Matematica Iberoamericana, 9:51-137, 1993; [12] J. Kovacevic and M. Vetterli. Nonseparable 
Multidimensional Perfect Reconstruction Filter-banks. IEEE Transactions on Information Theory, 
38:533-555, 1992; [11] W. He and M.J. Lai. Examples of Bivariate Nonseparable Compactly 
Supported Orthonormal Continuous Wavelets. In M. Unser, A. Aldroubi, A. Laine editor. Wavelet 
Applications in Signal and Image Processing IV, volume 3169 of Proceedings SPIE, pages 303-314, 
1997; [8] K. Grochenig and W. Madych. Multiresolution Analysis, Haar Bases and Self-Similar 
Tilings. IEEE Transactions on Information Theory, 38:558-568, 1992; [4] A. Ayache, E. Belogay, 
and Y. Wang. Orthogonal Lifting: Constructing New (Symmetric) Orthogonal Scaling Functions. 
2002; [5] E. Belogay and Y. Wang. Arbitrarily Smooth Orthogonal Nonseparable Wavelets in r^. 
SIAM Journal of Mathematical Analysis, 30:678-697, 1999; [1] E.H. Adelson, E. Simoncelli, and R. 
Hingoranp. Orthogonal Pyramid Transforms for Image Coding. In Visual Communications and Image 
Processing II, Volume 845 of Proceedings SPIE, pages 50-58, 1987; [18] E.P. Simoncelli, W.T. 
Freeman, E.H. Adelson, and J. P. Hager. Shiftable Multi-Scale Transforms. IEEE Transactions 
Information Theory, 38(2):587-607, 1992; and [19] J. Starck, E. J. Candes, and D.L. Donoho. The 
Curvelet Transform for Image Denoising. IEEE Transactions Image Processing, 1 1(6):670-684, 
2002), all dealing with different construction for performing multi-resolution analysis. Of note are 
the references of quincunx filtering schemes for performing multi-resolution analysis, which are 
stated in the specification as having "some continuity properties with respect to dyadic dilations or 
dilations induced by the Quincunx matrix only have been constructed in the past. All of them have 
no axial symmetries and are not smooth, except those constructed in E. Belogay and Y. Wang. 
Arbitrarily Smooth Orthogonal Nonseparable Wavelets in r^. SIAM Journal of Mathematical 
Analysis, 30:678-697, 1999, which can be made arbitrarily smooth, but are highly asymmetric." 

[20| M. Vetterli and J. Kovacevic. Wavelets and Subband Coding. Prentice Hall PTR, 
Englewood Cliffs, NJ. 1995); [13| D. Marr. Vision. A Computational hirestioatioii into the Human 
Representation and Pmccssino of Vi.sual Information. W. II. Freeman and Co., New York, N.Y., 
1982; and [14] D. Marr and F,. Ilidreth. The Theory of Fdge Detection. Proc. R. Soc. London B, 
207:187-217, 1980, deals with the fact that the human visual system critically depends on edge 
detection and with a specific construct to formula Laplacian filters to mimic the human visual 
system. 

[9| D. Han, D.R. Larson, M. Papadakis. and T. Stavropoulos. Multiresolution Analysis of 
Abstract Hilbert Spaces and Wandering Subspaces. In D. R. Larson L. Baggett, editor, r/zeFH/zmonrt/ 
and Harmonic Analysis of Wavelets and Frames, volume 247 of Com. Math., pages 259-2M.Amer. 
Math. Soc, 1999, relates to a proof that a is an injective homomorphism and a(G) is a subgroup 
ofG. 

[11] M. Pinsky. An Introduction to Fourier Analysis and Wavelets. 2001, relates to the proof 
of equation (7), which can be found in Lemma 2.5.1 of Pinsky. 

[17| M. Pinsky. An Introduction to Fourier Analysis and Wavelets. 2001 ; |3] G.E. Andrews, 
R. Askey, and R Roy. Special Functions. Number 71 in Encyclopedia of Mathematics. Cambridge 
University Press, 2000; and [21| G.N. Watson. A Treatise on the Theory of Bessel Functions. 
Cambridge Mathematical Library. Cambridge University Press, 1944, relate to details regarding 
Bessel functions and an extensive treatment of their main properties. 

[2] A. Aldroubi. Portraits of Frames. Proceeding of the American Mathematical Society, 
123: 1661-1668, 1995 and [10| D. Han and D.R. Larson. Frames, Bases and Croup Representations, 
volume 147 of Memoirs. American Mathematical Society, 2000, relate to supporting the statement 
in the specification that an arbitrary orthogonal projection R defined on a Hilbert space H maps 
every orthonormal basis of H onto a Parseval frame for R(H). 

The text in bold and italics is derived directly from the specification, where each reference was cited 
and can be foiond in the ft)llowing paragraphs: [0031], [0032], [0033], [0041], [0062], [0063], and 
[0078]. 



Applicants believe that the two Papadakis references were already submitted to the Examiner, 



but these two references form the background from which the method and system of this invention 
was derived, but in no way disclosed or taught now to go about the construction, which is the 
purpose of the present invention. 

Applicants attach a rule 132 Declaration concerning the Cohen et al. reference, in which Dr. 
Papakadis declares that the inventors did not rely on the Cohen et al. reference to derive the present 
invention. 



Specification 

The Examiner contends as follows: 

4. The specification is objected to because pages 26-27 contain a list of publications cited by 

the specification. According to the MPEP § 609 A(l) "the list may not be incorporated into the 
specification but must be submitted in a separate paper." Therefore, Examiner suggests Applicants file 
a separate information disclosure statement following the requirements of37 CFR 1.98(b), which 
requires a list of all patents, publications, or other information submitted for consideration by the 
Office. The list is improper and thus, must be removed from the specification. 

Applicants' attorney has never been required to remove as list of references from an 
application. Most of the references are not directly relevant to claim inventions of this invention, 
but they are merely direct readers to additional information about certain mathematical functions, 
derivations and/or understandings and understanding concerning human visual system. As shown 
above in the information disclosure statement section, the references are characterized as set forth 
expressly in the specification. Applicants have moved each reference to the place where it was cited 
in the specification so that readers can be directed to the articles or books that deal directly with the 
derivation or understanding being set forth at that location in the specification. 



Claim Rejections - 35 USC§ 101 
5. Claims 1-6, 13-24 stand rejected under 35 U.S.C. 101 as not falling within one of the four 
statutory categories of invention. 

The Examiner contends as follows: 

Federal Circuit precedent {In re Bilski, 88 USPQ2d 1385 (Fed. Cir. 2008)) requires that a 
statutory "process" under 35 U.S.C. 101 must (1) be tied to another statutory category (such as a 
particular apparatus), or (2) transform underlying subject matter (such as an article or material) to a 
different state or thing. While the instant claim( s) recite a series of steps or acts to be performed, the 
claim(s) neither transform underlying subject matter nor positively tie to another statutory category 
that accomplishes the claimed method steps, and therefore do not qualify as a statutory process under 
35 U.S.C. 101 (Claim 15' s recitation of a computer does not make the claim statutory because the 
computer is recited in the preamble and therefore, not given patentable weight.). 



Applicant's amendment to claims 1 and 4, adding the limitation "implemented on a computer" 
fails to overcome the 101 rejection above because the computer is recited in the preamble and thus, 
is not given any patentable weight. 

Applicants have amended the claims to add that the components are "in the form of software 
code executable on the computer," which is supported in the specification at least at paragraphs 12 
and 27. Applicants also state that the entire specification is directed to the construction of a system 
of software implemented on a computer designed to transform signals, information, data or images 
into wavelet decomposed and resolved sub-bands or resolution levels that can then be used to 
reconstruct a filtered resolution of the signals, information, data or images, which have improved 
properties and enhanced feature recognition. While the constructs that accomplish the signal, 
information, data or image transformation are computerized software derived from mathematical 
expression, the raw signals, information, data or images are not merely translated, scaled, rotated, 
or moved, but are filtered, decomposed and ultimately reconstructed into transformed signals, 
information, data or images with improved properties that greatly facilitate analysis. Filtering by 
definition results in a change in data. The original data is subjected to a set of filtering operation to 
enhance analysis of the data. In this case, the filtering is performed by completely isotropic, non- 
separable wavelets that are constructed from completely isotropic, non-separable windows, 
operations and filters, where the method can be extended to any dimension. Nothing in the prior art 
discloses or even suggests such a system or method. 

6. Claims 11-12 are rejected under 35 U.S.C. 101 because the claims merely recite 
mathematical functions. 

The Examiner contends as follows: 

Although a "computer" is recited in the claim, the recitation of the computer is in the 
preamble and thus, is not given patentable weight. 

Applicants have amended claims 1 1 and 12 to depend fi-om claim 1 and have amended claim 
1 to make it clear that the steps result in software code executable on the computer. Applicants, 
therefore, respectfiilly requests withdrawal of this rejection. 

7. Claims 13-14 are fitrther rejected under 35 U.S.C. 101. 

The Examiner contends as follows: 

The claims are nothing more than a mathematical algorithm. See MPEP 2106.02. 



Applicants have amended claims 13 and 14 to depend from claim 1 and have amended claim 
1 to make it clear that the steps result in software code executable on the computer. Applicants, 
therefore, respectfiilly requests withdrawal of this rejection. 



Claim Rejections - 35 USC§ 102 
8. Claims 7-10 stand rejected under 35 U.S.C. 102(b) as being anticipated by Bouchard et al., 
U.S. Patent No. 5,898,798 ("Bouchard"). 
The Examiner contends as follows: 

Referring to claim 7, Bouchard discloses a system for processing signals implemented on a 
computer comprising: 

a processing unit having encoded thereon a completely isotropic, non-separable 
ideal filter for frame multi-resolution analysis software including [col. 4, II. 13-35 
and col. 5, II. 57-67]: 

wavelets adapted to resolve a multidimensional signal into various resolution levels 
[col. 3, II. 29-39 and col. 4, II. 13-35], where the wavelets are derived from: 
isotropic, non-separable ideal windows or filters in a dimension 
greater than 1 , isotropic, non-separable low pass filters, isotropic, 
non-separable high pass filters and isotropic, non-separable 
filters that cover a desired frequency range or plurality of 
frequency ranges; and isotropic, non-separable frame scaling 
functions where translations of the frame scaling functions form 
a frame [col. 4, II. 13-35. Note the bidimensional non-separable 
low and high-pass filters and scaling functions that allow the 
image to be processed in an isotropic way.]; 
where the system resolves or decomposes multidimensional signals, data, information, or images into 
a plurality of non-overlapping sub-bands sets or resolution levels with the at least one isotropic, non- 
separable wavelet improving analysis efficiency and improving analysis of more complex 
multidimensional signals, data, information or images [col. 3, II. 29-39 and col. 4, II. 13- 35]. 

Quincunx filters are non-separable, but are not isotropic. In fact, Bouchard entire system is 
built arovind preparing quincunx filter banks that are bidimensional, but are not isotropic as is clearly 
shown in Figijre 3 and the associated text of Bouchard. Bouchard's bidimensional non-separable low 
and high-pass filters are bases on a quincunx grid and are rotational invariant for rotation angles of 
pi/4. Bouchard at Col. 4, 13-35. Bouchard designed the bidimensional non-separable low and high- 
pass filters to accommodate a specific rotational symmetry. While the bidimensional non-separable 
quincimx wavelets have limited rotational symmetry that is designed into the wavelets when they 
are formed, the bidimensional non-separable quincunx wavelets are not isotropic. While it is 
true that quincunx wavelets can be designed to be nearly isotropic, they are not isotropic (see M. 
Feilner, D.V. De Ville, and M. Unser, "An Orthogonal Family of Quincunx Wavelets With 



Continuously Adjustable Order," IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 14, 



NO. 4, APRIL 2005 attached). Applicants also directly and expressly address quincunx wavelets 
in specification stating: 

[0032] In this invention, we construct non-separable Sliannon-like GFMRAs of i^(R") 
wliose scaling functions are radial and are defined with respect to certain unitary systems, which we 
will later introduce. We also derive certain of their associated frame muhiwavelet sets. Our 
construction is the first of its kind. Scaling functions that are radial have not been constructed in the 
past. However, certain classes of non separable scaling functions in two dimensions, with some 
continuity properties with respect to dyadic dilations or dilations induced by the Quincunx 
matrix only have been constructed in the past (e.g., [7, 12, 11, 8], |4]). All of them have no axial 
symmetries and are not smooth, except those constructed in |5], which can be made arbitrarily 
smooth, but are highly asymmetric. Another construction in the spirit of digital filter design, but not 
directly related to wavelets can be found in [1] and [18], The latter construction and this of curvelets 
(e.g., see [19]) share two properties of our Radial GFMRAs: the separability of the designed filters 
with respect to polar coordinates and the redundancy of the induced representations. However, our 
construction in contrast to those due to Simoncelli et. al. and to Starck et al. are in the spirit of classical 
multiresolution analysis and can be carried out to any number of dimensions and with respect to a 
variety of dilation matrices. 

Specification at Paragraph [0032] (emphasis added). 

The present invention relates to wavelets that are isotropic, not nearly so. 

Moreover, quincunx wavelets and quincunx multi resolution analysis is limited to 2 
dimensions only and is not extendable to 3 or more dimension, save for accumulating a set of 2D 
slices of a 3D image, where each 2D slice is analyzed using a quincunx formalism (see D.V. De 
Ville, Thierry BIu, and Michael Unscr, "On the Multidimensional Extension of the Quincunx 
Subsampling Matrix," IEEE SIGNAL PROCESSING LETTERS, VOL. 12, NO. 2, FEBRUARY 
2005 attached). The present invention utilized completely isotropic, non-separable windows, 
operators, filters and wavelets to perform multiresolution analysis without any dimensional 
constrains. 

Because Bouchard does not disclose or teach completely isotropic, non-separable windows, 
operators, filter and wavelets of this invention, Bouchard cannot anticipate claim 7 as claim 7 relates 
specifically to a system based on completely isotropic, non-separable windows, operators, filter and 
wavelets. Moreover, Bouchard does not even suggestion or directly an ordinary artisan to fiinction 
that are completely isotropic and non-separable as Bouchard is tied to quincunx bidimensional 
wavelets. Applicants, therefore, respectfiiUy request withdrawal of this rejection. 



The Examiner contends as follows: 

Referring to claims 8-10, Bouchard further discloses that each isotropic, non-separable high pass and 
each isotropic, non-separable low pass filter comprise: a plurality of isotropic, nonseparable high pass 
and isotropic, non-separable low pass components, each component including at least one isotropic. 



non-separable relative low pass subcomponent and at least one isotropic, non-separable relative high 
pass subcomponent [col. 4, II. 13-35 and fig. 3]. 

Applicants reassert their argument concerning Bouchard here. Bouchard does not disclose 
isotropic, non-separable high pass and low pass filters, especially filters constructed from isofropic 
components. The Bouchard filters are rotational invariant with respect to a specific angle meaning 
that they will have a directional bias and the Bouchard filters are imitated to 2D analj^is. The 
present invention utilizes completely isotropic windows, operators, filters and wavelets to 
decompose or resolve signal or image into subbands without regard to the nature of the image as the 
windows, operators, filters and wavelets are all isofropic. The present system then reconstructs an 
filtered and enhanced image from the decomposition. 

Because Bouchard does not disclose or teach completely isotropic, non-separable windows, 
operators, filter and wavelets of this invention, Bouchard cannot anticipate claim 7 as claim 7 relates 
specifically to a system based on completely isotropic, non-separable windows, operators, filter and 
wavelets. Moreover, Bouchard does not even suggestion or directly an ordinary artisan to function 
that are completely isotropic and non-separable as Bouchard is tied to quincunx bidimensional 
wavelets. Applicants, therefore, respectfiiUy request withdrawal of this rejection. 

Having fiilly responded to the Examiner's Non-Final Office Action, Applicants respectfully 
virge that this application be passed onto allowance. 

If it would be of assistance in resolving any issues in this application, the Examiner is kindly 
invited to contact applicant's attorney Robert W. Strozier at 713.977.7000 

The Commissioner is authorized to charge or credit Deposit Account 501518 for any 
additional fees or overpayments. 

Respectfully submitted. 

Date: 8 September 2009 

/Robert W. Strozier/ 

Robert W. Strozier, Reg. No. 34,024 



An Orthogonal Family of Quincunx Wavelets 
With Continuously Adjustable Order 

Manuela Feilner, Dimitri Van De Ville, Member, IEEE, and Michael Unser, Fellow, IEEE 



Abstract — We present a new family of two-dimensional and 
three-dimensional orthogonal wavelets which uses quincunx sam- 
pling. The orthogonal refinement filters have a simple analytical 
expression in the Fourier domain as a function of the order A, 
which may be noninteger. We can also prove that they yield wavelet 
bases of L^iU^) for any A > 0. The wavelets are fractional in 
the sense that the approximation error at a given scale a decays 
like 0(a*); they also essentially behave like fractional derivative 
operators. To make our construction practical, we propose a fast 
Fourier transform-based implementation that turns out to be 
surprisingly fast. In fact, our method is almost as efficient as the 
standard Mallat algorithm for separable wavelets. 

Index Terms — McClellan transform, nonseparable filter design, 
quincunx sampUng, wavelet transform. 



Tllli CiRliAT majorily of wavelet bases that are currently 
used Cor image processing are separable. There are two pri- 
mary reasons for this. The first is convenience, because wavelet 
theory is most developed in one dimension and that these results 
are directly transposable to higher dimensions through the use 
of tensor product basis functions. The second is efficiency be- 
cause a separable transform can be implemented by successive 
one-dimensional (1-13) processing of the rows and columns of 
the image. The downside, however, is that separable transforms 
tend to privilege the vertical and horizontal directions. They also 
produce a so-called "diagonal" wavelet component, which does 
not have a straightforward directional interpretation. 

Nonseparable wavelets, by contrast, offer more freedom and 
can be better tuned to the characteristics of images [ 1 ] , [2] . Their 
less attractive side is that they require more computations. The 
quincunx wavelets are especially interesting because they can 
be designed to be nearly isotropic [3]. In contrast with the sep- 
arable case, there is a single wavelet and the scale reduction 
is more progressive: a factor \/2 instead of 2. The preferred 
technique for designing quincunx wavelets with good isotropy 
properties is to use the McClellan transform [4] to map 1-D 
biorthogonal designs to the multidimensional case. Since this 
approach requires the filters to be symmetric, it has only been 
applied to the biorthogonal case because of the strong incen- 
tive to produce filters that are compactly supported [5]-[8]. One 
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noteworthy exception is the work of Nicolier et al. who used 
the McClellan transform to produce a quincunx version of the 
Battle-Lemarie wavelet filters [9]. However, we believe that 
their filters were truncated because they used a representation 
in terms of Tchebycheff polynomials. 

In this paper, we construct a new family of quincunx wavelets 
that are orthogonal and have a fractional order of approxima- 
tion. The idea of fractional orders was introduced recently in 
the context of spline wavelets for extending the family to non- 
integer degrees [10]. The main advantage of having a continu- 
ously varying order parameter — ^not just integer steps as in the 
traditional wavelet families — ^is flexibility. It allows for a con- 
tinuous adjustment of the key parameters of the transform, e.g., 
regularity and localization of the basis functions. The price that 
we are paying for these new features — orthogonality with sym- 
metry as well as fractional orders — ^is that the filters can no 
longer be compactly supported. We will make up for this hand- 
icap by proposing a fast Fourier transform (FFT)-based imple- 
mentation which is almost as efficient as Mallat' s algorithm for 
separable wavelets [11]. 

II. Quincunx Sampling and Filterbanks 

First, we recall some basic results on quincunx sam- 
phng and perfect reconstruction filterbanks [12]. The quin- 
cunx sampling lattice is shown in Fig. 1. Let x[k] with 
k = (ki,k2) 6 1^ denote the discrete signal on the ini- 
tial grid. The two-dimensional (2-D) z-transform of x[k{\ is 
denoted by = Xl^g^^ •■^t^]-^^' where 5* = z^'^Z2^. 

The continuous 2-D Fourier transform is then given by 
^(e^'") = Eftez^ with a; = (a;i,a;2), and, 

finally, the discrete 2-D Fourier transform for given on 
an X TV grid {ki,k2 = 0, 1, . . . , AT - 1) by X[n\ = 
Efeez2 a:[^]e-J'2-<"'^>/^', with m, na = 0, 1, . . . , - 1. 

Now, we write the quincunx sampled version of x[k] as 



[x]i-D[k] = x[Dk], where 



-0 -0 



(1) 



Our down-sampling matrix D is such that = 21. The 
Fourier-domain version of (1) is 

WlD[fc] ^\[X (e^-°-"^) + X (2) 

where tt = (7r, tt). 

The upsamphng is defined by 



k]= I x[T)-^k], when ki + k2 h 
0, elsewhere 
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(a) (b) 

Fig. 1. (a) Quincunx lattice and (b) the corresponding Nyquist area in the frequency domain. 
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H{z) 
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-Xi[k] 



and its effect in the Fourier domain is as follows: 

MTD[fc] — (4) 

If we now chain the down-sampling and up-sampling operators, 
we get 



Mj,DTD[fc] = 



r H{m{^) + G{z)G{z) = 2 
\ H{-z)H{z) + G{-z)G{z) = 0 



(7) 



where H and G (respectively H and G) are the transfer func- 
tions of the synthesis (respectively analysis) filters. In the or- 
thogonal case, the analysis and synthesis filters are identical up 
to a central symmetry; the wavelet filter G is simply a modu- 
lated version of the low-pass filter H. 

III. Fractional Quincunx Filters 
To generate quincunx filters, we will use the standard ap- 
proach which is to apply the diamond McClellan transform to 
map a 1-D design onto the quincunx structure. 



filterbank on a quincunx lattice. 



A. New 1-D Wavelet Family 



As starting point for our construction, we introduce a new 1 -D 
family of orthogonal filters 

V2{z + 2 + z-')^ 



^(z + 2 + z-^)^ + i~z + 
v/2(2-F2cosa>)7 



(6) 



Since quincunx sampling reduces the number of image samples 
by a factor of two, the corresponding reconstruction filterbank 
has two channels (cf. Fig. 2). The low-pass filter H reduces the 
resolution by a factor of \/2; the wavelet coefficients correspond 
to the output of the high-pass filter G. 

Applying the relation (6) to the block diagram in Fig. 2, it is 
easy to derive the conditions for a perfect reconstruction 



which is indexed by the continuously-varying order parameter 

These filters are symmetric and are designed to have zeros 
of order A at ^ = —1; the numerator is a fractional power 
of -I- 2 -I- ^;~^) (the simplest symmetric refinement filter of 
order 2) and the denominator is the appropriate /2-orthonormal- 
ization factor. By varying A, we can adjust the frequency re- 
sponse as shown in Fig. 3. As A increases, H\{z) converges 
to the ideal half-band low-pass filter. Also note that these fil- 
ters are maximally flat at the origin; they essentially behave like 
Hx{u))/-V2 = 1 + O(w^) as cj ^ 0. Their frequency response 
is similar to the Daubechies' filters with two important differ- 
ences: 1) the filters are symmetric and 2) the order is not re- 
stricted to integer values. 

We can prove mathematically that these filters will generate 
valid 1 -D fractional wavelet bases of L2 similar to the fractional 
splines presented in [ 10] . The order property (here fractional) is 
essential because it determines the rate of decay of the approx- 
imation error as a function of the scale. It also conditions the 
behavior of the corresponding wavelet 0 which will act like a 
fractional derivative of order A. In other words, it will kill all 
polynomials of degree n < [A — 1] ; i.e. 



x"i,x{x)da; = 0. 



(9) 
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Fig. 3. Frequency responses of the orthogonal refinement filters for 
A = 1,...,100. 

B. Corresponding 2-D Wavelet Family 

Applying the diamond McClellan transform to the filter 
above is straightforward; it amounts to replacing cosw by 
(1/2) (cos wi + costi;2) in (8). Thus, our quincunx refinement 
filter is given by 

jj ^g?^^ _ %/2(2+cosc<;i +008 0)2)^ 

\/ (2+cosa;i+cosa;2)'^ + (2 — coswi — cos 0^2)'^ 
(10) 

This filter is guaranteed to be orthogonal because the 
McClellan transform has the property of preserving biorthog- 
onality. Also, by construction, the Ath order zero at a; = tt 
gets mapped into a corresponding zero at cD = (tTjTt); 
this is precisely the condition that is required to get a 2-D 
wavelet transform of order A. Also, note the isotropic be- 
havior and the flatness of H\{e^'^) around the origin; i.e., 
Hxie^^)/V2 = l + 0(||d;||^)forw ^ O.Fig.4showscontour 
plots of the scaling filter for several choices of the order A. 

The orthogonal wavelet filter is obtained by modulation 

Gx(z) = z,Hx(-z-^). (11) 

The corresponding orthogonal scaling function (p\{x) is de- 
fined implicitly as the solution of the quincunx two-scale rela- 
tion 

M^v) = V2J2 hx[k]ipxiDx-k). (12) 

Since the refinement filter is orthogonal with respect to the quin- 
cunx lattice, it follows that (fi\{x) G ^2(1^^) and that it is or- 
thogonal to its integer translates. Moreover, for A > 0, it will 
satisfy the partition of unity condition, which comes as a direct 
consequence of the vanishing of the filter at (wi, W2) = (tt, tt). 
Thus, we have the guarantee that our scheme will yield orthog- 
onal wavelet bases of L2{U'^). The underlying orthogonal quin- 
cunx wavelet is simply 

Va(^) = V2Y,ax [k]>Px (D.T - k) . (13) 
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Fig. 4. Contour plots of the low-pass filters Hx{e'") for various values of the 
order parameter A. (a) A = 1. (b) A = y/2. (c) A = tt. (d) A = 10. 

Given the behavior of Hx{pJ^) at cJ = 0, we also have V-'a(w) oc 
||a;||'^, and, as such, the wavelet behaves as the Ath order differ- 
entiator for low frequencies [13]. The vanishing moment prop- 
erty in the 2-D case becomes 

j x1'x2'tpx{x)dx, for ni+n2< (14) 

Fig. 5 shows the wavelet 'ipx{x) for various choices of the order 
A. Note that the wavelet is centered around (1/2, 1/2). As illus- 
trated by these plots, the wavelets clearly gets smoother as A 
increases. However, a mathematical rigorous estimation of their 
regularity is beyond the scope of this paper. 

IV. Implementation in Fourier Domain 

The major objection that can be made to our construction is 
that the filters are not FIR and that it may be difficult and costly 
to implement the transform in practice. We will see here that we 
can turn the situation around and obtain a very simple and effi- 
cient algorithm that is based on the FFT, following the idea of 
[14]. Working in the frequency domain is also very convenient 
because of the way in which we have specified our filters [see 
(10) and (11)]. Implementations of the wavelet transform for the 
quincunx subsampling matrix using FFTs have been proposed 
before [9], [15]; our algorithm is another variation, which min- 
imizes the number and size of FFTs and seems to be faster. 

First, let us assume that the image size is TV x TV. Now, we will 
describe the decomposition part of our algorithm which corre- 
sponds to the block diagram presented in Fig. 6, where we have 
pooled together two levels of the decomposition. The initializa- 
tion step is to evaluate the FFT of the initial input image x[k] and 
to precompute the corresponding sampled frequency responses 
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(a) X=\/2 (b) A - TT (c) A = 10 

Fig. 5. Surface plots of the wavelets V'a for various values of the order parameter, (a) A = \/2. (b) A = tt. (c) A = 10. 



Xi[k]- 



GA^ \ -{i}^Vi+2[k] 



y'i+i[k] 



21 

-<i>- 



Fig. 6. Analysis part of the 2-D quincunx wavelet transform for two 



of the analysis filters H[7f] and G[n] using (10) and (11). We also 
precompute the rotated version of the filters, denoted as Hp[n] 
and Gp[n], that can be obtained as 

Hp[n] = H [Dn mod {N, N)] (15) 

G„[n] =G[D»niod(Af,AO]. (16) 

Let us now consider the 2-D FFT of the input, given by 



Xi[n] = '^x,[k]c-^^^^^, forni,n2 = 0, 1, . . . , - 1. 

(17) 

Globally, at the end of the process, tiie output variables are 
the quincunx wavelet coefficients yi[k], 2/2 [fc], • • • ^yj[k], and 
e.g., as shown in Fig. 7(a). Their Fourier transforms for 
the odd iterations are derived from the auxiliary N x N signals 
(see also Fig. 6) 



fc 



Down and up sampfing with D in the first iteration step intro- 
duces zeros in the space domain while it preserves the size of 
Y,-_^_-^ [n] . However, it implies some symmetiy /redundancy in fre- 
quency domain. Therefore, only half of the coefficients needs 



to be computed which saves operations. The reduced signal 
Y^_^-^ [fc] and its corresponding low-pass signal are obtained by 



N 

I 2 ' 2 



)]) 



Xl+An'] = \ (H[n']XAn'] + H [rf + |) 



(18) 
(19) 



where n' G [0, {N/2) - 1] x [0, TV - 1]. 

To generate the signal of (19) in the way that is de- 

picted in Fig. 7(a) with every second row shifted by one pixel, 
we separate the image in even ('(/t+i,even) and odd ('(/i+i,odd) 
rows already in the Fourier domain, using the auxiUary variable 
Z[m] 

Z[m] = Y/^,[m] + y/+i [m + (o, y)] 



Zi+l[k] =!/i+l,eve„[fc] + jVi^ 



^i+l,odd[fc] <^23) 

with rn G [0, {N/2) - 1]^. The sum in the real part {Y/^^ [rn] + 
Y/^-^[m + (0,iV/2)]) represents downsampling by two in the 
vertical direction, keeping aU the even rows, whereas the sum 
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(a) (b) 

Fig. 7. Wavelet coefficients for the quincunx subsampling scheme can be arranged in two ways. An example for 7 = 4 iterations, (a) Compact representation. 

(b) Classic representation. 

in the imaginary part represents tlie odd rows. In the space do- 
main, we alternate the rows 2A;2 + 1] = Ro{2;[fc]} and 
2k<2\ = Im{2;[A;]}. Since z[k] is four times smaller than 
y'i+i [^]' we save computations with the reduced-size IFFT. 

Instead of rotating the frequency variables after each itera- 
tion, we use the precomputed rotated version of the filters (i.e.. 
Hp and Gp), which we apply at all even iterations. In this way, 
we also save two rotations per iteration in the frequency domain. 

The Fourier transforms of the output for the even iterations 
are 



position algorithm using up sampling, instead. For instance, the 
synthesis counterpart of (25) and (26) is 

' ''■"2 + y ^ j = ^■'■+2 ' ''■"a] , 
^/+2 |'mi,m2+ ^y^l =Yi+2[mi,m2\, 

+ ^:+2[™i>"2]Gp[mi,n2]. 



y,+2[m] = 5:?y.+2[fc]e"''^ 
k 

N 

formi,m2 = 0, - 1. (24) 

They are computed by 

^i+2[m] = i (^Hp[m]Xl+,[7n] + Hp + (^0, y)] 

Yi+^[fh] = \ {Gp[m]X[^Am] + Gp \m. + [q., 

xX^+i [m+ (|o,y)]). (26) 

The process is then iterated untU one reaches the final resolution. 
When the last iteration is even, we lower the computation costs 
with the FFT by utilizing its imaginary part 

z[k]=Y, {Xi+2 [m] + jFi+s H) e V' ' (27) 



where .Ti+2[fc] = Re{z[fc]} andyi+2[A^] = Im{z[fc]}. 

Obviously, as the resolution gets coarser after each iteration, 
the Fourier transforms of the filters need not be recalculated; 
they are simply obtained by down-sampling the previous arrays. 

The synthesis algorithm operates according to the same prin- 
ciples and corresponds to the flow graph transpose of the decom- 



V. I EXPERIMENTS 
A. Benchmark and Tesliiig 

We have implemented two versions of the algorithm, based 
on Java and Matlab. For the Matlab version, we report compu- 
tation times below 0.8 s for 16 quincunx iterations of a 256 x 
256 image on an Apple G4 700 MHz desktop; the decompo- 
sition is essentially perfect with a reconstruction error below 
10"-^^ RMS. The method is generic and works for any set of 
filters that can be specified in the frequency domain, indepen- 
dent of their spatial support (or infinite spatial support, such 
as in our case). As a comparison, the Matlab implementation 
available in the latest Wavelet Toolbox [16] for the Daubechies 
9/7 filters (used in JPEG 2000) applied to the same image and 
for an equivalent of eight separable iterations, takes about 1 .7 s. 
For N datapoints, the complexity of our approach boils down to 
0{N log N) for theFFT-based implementation, versus 0{NB) 
for the spatial-domain implementation, where B is related to the 
filter support. The exact tradeoff will depend on the image size 
and the filter size. However, taking into account the benchmark 
measures and its flexibility, we beUeve that the FFT-based im- 
plementation deserves consideration for a broad class of appli- 
cations. 

We also provide an applet written in Java, which makes 
it possible to run the algorithm over the Internet, at the site 
http://bigwww.epfl.ch/demo/jquincunx/. A screen shot of this 
applet is presented in Fig. 8. 

Two examples of fractional quincunx wavelet decomposi- 
tions with A = \/2 and A = tt are shown in Fig. 9. Note how 
the residual image details are more visible for the lower value 
of A. The larger A reduces the energy of the wavelet coefficient, 
but this also comes at the expense of some ringing. Thus, it is 
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Fig. 8. Applet of the Fourier-based implementation of the q. 



'ailable on the site: http://bigwww.epfl.ch/demo/jquincunx/. 




(a) input image 



IX wavelet transforms with four 



/2 (c) A - 7r 

(a) Original test image, (b) A = \/2 • (c) A = tt 



convenient to have an adjustable parameter to search for the best 
tradeoff. 

An advantage of the present approach is that the filters for 
small A are nearly isotropic; this is the reason why the wavelet 
details in Fig. 9 do not present any preferential orientation. The 
degree of isotropy of the various lowpass filters can be seen 
from Fig. 4. The shape of the contour plots of the low-pass filter 
Hx {e->'^ ) confirms that the degree of isotropy is the best for small 
values of A. At the other extreme, when A — *■ oc, H\{e^'^) tends 
to the diamond-shaped ideal filter. 

Another nice feature of the algorithm is that the computa- 
tional cost remains the same irrespective of the value of A. 

B. Dependence of the Order Parameter 

The usefulness of a tunable order parameter is demonstrated 
in the following experiment: we apply the quincunx transform 
to the test image "cameraman" [see Fig. 10(a)] and reconstruct 
using only 15% of the largest coefficients. Then the SNR is 
measured depending on the order parameter. The plot in Fig. 1 1 
shows how the SNR changes according to the order A; the 
optimum, indicated by the circle, is achieved for A = 2.5. 
Fig. 10(b) and (c) shows the reconstructions for the optimal 
order and an order too high. The last one gets penalized by the 
introduction of ringing artefacts around the edges. We also plot 



the SNR curves for 20% and 25% of the coefficients. The same 
type of qualitative behavior holds for other images. 

C. Approximation Properties 

The main differences between the quincunx and the conven- 
tional separable algorithm is the finer scale progression and the 
nonseparability. To test the impact that this may have on com- 
pression capability, we compared the approximation qualities of 
both approaches. Since the wavelet transform is orthogonal, the 
approximation error (distortion) is equal to = ||x — = 
1 1 J/ — y 1 1 ^ , where y are the wavelet coefficients of the input image 
x; X is the reconstructed image obtained from the quantized — or 
truncated — wavelet coefficients y. Also, in the space domain 
is equivalent to the sum of squares of discarded wavelet coeffi- 
cients [17]. 

1 ) Linear Approximation: In classical rate-distortion theory, 
the coefficients are grouped into channels and coded indepen- 
dently. In the orthogonal case, is equivalent to the differ- 
ence between the signal's energy and the energv of the recon- 
structed signal - = - ||;rf = - (\\:rj\\^ + 

J2'j=i liyjlP)- The distortion across N channels with variance 



D = N -C- 2-2« . 
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order parameter X 



Fig. 11. Relation between the order parameter A and the SNR of the 
reconstructed image (test image "cameraman") using only the largest 
coefficients. The full hne, dashed hne, and dotted line correspond, respectively, 
to 25%, 20%, and 15% of the largest coefficients. 

where C is a constant, R is the mean rate, and p is the geometric 
mean of the subband variances 

p=(n-^') • (29) 

When p is small, the distortion is small, as well. What this means 
quaUtatively is that the wavelet transform which has the larger 
spread in the variances will achieve the better coding gain [12]. 
The linear approximation subband coding gain for sample-by- 
sample quantization (PCM) is described by 

To better illustrate this issue, we have decomposed the test 
image "cameraman" for the maximal number of iterations, both 
for the quincunx and the separable case as shown in Fig. 12. 
The order was fixed (i.e., A = 4) for our method and for the or- 
thogonal separable approach (corresponding to the commonly 



used degree parameter a = 3 for the underlying B-splines). In 
Fig. 13(a), we compare the energy packing properties of both 
decompositions for linear approximation. "Fnergy packing" 
refers to the property that the more the first coefficients contain 
energy, the better the DWT yields compression. We start to 
sum up the energy of the subbands with the lowest resolution. 
Each step of the stairs represents a subband. ' The first subbands 
of the quincunx decomposition report higher energy packing 
than the separable case, but the overall coding gain is slightly 
better for the separable case than the quincunx case (47.69 
versus 45.23). Fig. 13(c) shows similar results for the "Lena" 
test image. 

Since the branches are orthogonal, the transformation that 
provides the maximum energy compaction in the low-pass 
channel is also the one that results in the minimum approxi- 
mation error [17]. Since most images have a power spectrum 
that is roughly rotationaUy invariant and decreases with higher 
frequencies, separable systems are usually not best suited for 
isolating a low-pass channel containing most energy and having 
high-pass channels with low energy. In contrast, a quincunx 
low-pass filter will retain more of the original signal's energy 
[12]. 

Consequently, the type of images that benefit the most from 
the quincunx scheme have a more isotropic spectrum. For 
example, for the well-known zoneplate test image of Fig. 14(a), 
the coding gain of quincunx scheme is about 20% better than 
the one obtained by the separable scheme (4.30 versus 3.64). 
Also, the quincunx scheme gives better energy compaction for 
textures of highly isotropic nature (and as such a higher coding 
gain). Two such examples of the Brodatz textures are shown 
in Fig. 14(b) and (c), corresponding to a coding gain of 13.67 
versus 12.45 and 12.04 versus 9.62, respectively. On the other 
hand, a separable treatment leads to a better energy compaction 
for the texture shown in Fig. 15 (8.78 versus 15.48). Other 
authors have also found that texture analysis using the quin- 
cunx scheme improves the results as compared to the separable 
scheme [18]. 

'A quincunx wavelet decomposition with J iterations generates N = J 
channels, while a separable wavelet decomposition with J iterations results into 
= 3/ -I- 1 channels. 
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number of coefficients 



number of coefficients 



Fig. 13. Comparison of energy-compaction property between ttie quincunx and ttie separable case of image decomposition (as stiown in Fig. 12). (a) and 
(c) Linear approximation depending on number of coefficients (in log, grouped per subband), respectively, for "cameraman " and "Lena." (b) and (d) Nonlinear 
approximation depending of the n largest coefficients (in log), respectively, for "cameraman" and "Lena." The quincunx scheme yields better results for a low 
number of coefficients. In the case of "Lena," the separable scheme performs better than the quincunx one over most of the range. 



2) Nonlinear Approximation: A more recent trend in 
wavelet theory is the study of nonlinear approximation. In this 
case we do not take the "n-first" coefficients, but the "n-Iargest" 
coefficients to approximate a signal with 71 coefficients. This 
yields better energy packing, since in the wavelet domain 
the "«-first" coefficients are not necessarily the largest one, 



especially along the position indices [19]. The distortion is 
described by [20] 



D^ = \\y-yK(T)\\' = E l2/Nr (31) 

\y\<T 
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Fig. 16(a). Let .7;[fc] denote the discrete signal on the initial grid. 

Then, its quincunx sampled version, following [6], is 



[.r] ID {^\ = ■r[Dk] ■ where D = 



Our down-sampling matrix D is such that D"* = 21 and 
I det D| =2. The Fourier-domain version of this formula is 
similar to the 2-D case 



re (Brodatz D58) that is better suited for a separable 



Moreover, it can be shown that 

D(xC-{N{T))-i 



(34) 



where tt = (tt, tt, tt). 

The implementation for the 3-D case goes as follows. The 
output variables are the discrete Fourier transforms of the 
wavelet coefficients 



(32) 



when the smoothness of y is measured by its inclusion in some 
critical Besov space B^{L'i{I)) with 1/q = (7/2) + (1/2), 
roughly speaking, when y is a fiinction with 7 derivatives in 
LHI) [20], [21]. 

For the nonlinear approximation, the quincunx scheme also 
yields a better approximation than the separable one for a small 
n in many cases. Fig. 13(b) represents the energy depending on 
the n largest coefficients (in log). 

VL Extension to Three Dimensions 

The extension of quincunx sampling to three dimensions 
is rather straightforward. First, the filters are obtained by re- 
placing cos w by (1/3) (cos uji + cos uj2 + cos wg) in (8). Next, 
the quincunx sampling lattice for three dimensions is shown in 



k 

for ni, 712, n3 = 0, 1, . . 
yi+2[n] = 5^t/,+2[fc]e--''^ 



The coefficients themselves are recovered by inverse FFT. The 
Fourier transforms after the first level of filtering are given by 

-4-(f'f^f)])- « 
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After the second level of filtering, we have 

xAVi.[n+(o.f...)j) (40) 




xX,+i ^n+(o, (41) 

Note that these are computed at the resolution of the input. The 
size reduction only takes place during the third step 

Xi+M] = \ (ffpp[m]X,:+2[rn] + H„ \m + (^^, 0, o)] 
xX.,.[r.+ (f,0,0)]) (42) 

xXi+2[m+(^^,0,o)]) (43) 

where Hp[m] = H[Y>rn mod {N,N,N)] and Hpp[rn] = 
H[T)'^r7i mod J^N.N.N)]. Analogously, we^have that: 
Gj,H = G[Yi77i mod {N,N,N)] and Gpp[m] = 
G[D'^m mod {N, N, N)]. 

Fig. 16(b) shows how the coefficients can be arranged in a 
nonredundant way inside the cube. Note that the size of the FFTs 
for the 3-D implementation can be further reduced by taking 
into account the subsampled arrangement in the Fourier domain. 
Again, the rotated filters Hp, Hpp, Gp, and Gpp are precom- 
puted. 



A. Approximation Properties in Three Dimensions 

We compared the compression capability for the quincunx 
and the separable scheme applied to 3-D data, similar to the type 
of experiments that are described for two dimensions in Sec- 
tion V-C . Fig. 1 7 shows the results for a spiral CT dataset of part 
of a human spine. The linear approximation quaUty is shown in 
Fig. 17(b). The separable scheme takes much advantage of the 
availibily of many small (i.e., seven for each iteration) bandpass 
subbands, as compared to the quincunx scheme. To illustrate 
this point, we have grouped the bandpass subbands for the sep- 
arable case together in one single bandpass in Fig. 17(c). For 
nonlinear approximation, both schemes perform similarly with 
a small advantage for the separable one, as shown m Fig. 17(d). 
If the dataset contains more (isotropic) high-frequency compo- 
nents, the breakpoint between the quincunx and the separable 
case shifts to the right. 

The main advantage of the 3-D quincunx scheme is in ap- 
pUcations that can benefit from the (much) slower scale pro- 
gression. One example is the statistical analysis of brain activity 
using functional magnetic resonance imaging (fMRI). Here we 
show an example using the classical wavelet-based approach for 
detecting activity, using the linear model analysis and the t test 
in the wavelet domain for a 3-D dataset (64 x 64 x 64) with 
an auditory stimulus [22]; we refer to [23] for more details. We 
compared the use of the 3-D dyadic separable wavelet decom- 
position based on orthogonal hnear B-sphne wavelets versus 
our 3-D quincunx wavelets (same order). The parameter maps 
where obtained using tiie same threshold after reconstruction 
(5% of the maximal parameter value). The number of detected 
voxels, and as such the sensitivity of the approach, is almost 
10% higher (578 versus 536) when we use the 3-D quincunx 
DWT, which confirms that the slower scale progression im- 
proves the quality of the results. Fig. 18 shows the detected ac- 
tivation patterns around the auditory cortex (slice 33). 

Other potential applications might include image analysis and 
3-D feature detection. 
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Fig. 17. (a) Slice of an spiral CT dataset of part of a human spine (courtesy and copyright of Ramani Pichumani, Stanford University School of Medicine), 
(b) Linear approximation, for the separable case each bandpass subband is considered independendy. (c) Linear ^proximation, for the separable case the bandpass 
subbands are grouped together into one single subband. (d) Nonlinear approximation. 



VII. Conclusion 
We have introduced a new family of orthogonal wavelet trans- 
forms for quincunx lattices. A key feature is the continuously 
varying order parameter A which can be used to adjust the band- 
pass characteristics as well as the locahzation of the basis func- 
tions. 

We have also demonstrated that these wavelet transforms 
could be computed quite efficiently in two and three dimensions 



using FFTs. This should help dispel the commonly held belief 
that nonseparable wavelet decompositions are computationally 
much more demanding than the separable ones. 

Because of their nice properties and their ease of implementa- 
tion, these wavelets present an alternative to the separable ones 
that are being used in a variety of image processing applications 
(image analysis, image enhancement, filtering and denoising, 
feature detection, texture analysis, and so on). 
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On the Multidimensional Extension of the Quincunx 
Subsampling Matrix 
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Abstract — The dilation matrix associated witli tlie tliree-dimen- 
sional (3-D) face-centered cubic (FCC) sublattice is often consid- 
ered to be the natural 3-D extension of the two-dimensional (2-D) 
quincunx dilation matrix. However, we demonstrate that both di- 
lation matrices are of different nature: while the 2-D quincunx ma- 
trix is a similarity transform, the 3-D FCC matrix is not. More 
generally, we show that is impossible to obtain a dilation matrix 
that is a similarity transform and performs downsampling of the 
Cartesian lattice by a factor of two in more than two dimensions. 
Furthermore, we observe that the popular 3-D FCC subsampling 
scheme alternates between three different lattices: Cartesian, FCC, 
and quincunx. The latter one provides a less isotropic sampUng 
density, a property that should be taken into account to properly 
orient 3-D data before processing using such a subsampling ma- 
trix. 

Index Terms — Nonseparable design, two-channel filterbanks, 
wavelet decomposition, 2-D quincunx sampling, 3-D FCC sam- 
pling. 



I. Introduction 

THE QUINCUNX sampling scheme is an interesting con- 
figuration that is commonly used for designing critically 
sampled filterbanks and discrete wavelet transforms in two di- 
mensions [l]-[7]. In contrast to the dyadic separable case, the 
quincunx dilation matrix leads to a two-channel filterbank archi- 
tecture — or equivalently one single wavelet — and reduces the 
scale more progressively: a factor \/2 instead of 2. We briefly re- 
view the two-dimensional (2-D) quincunx scheme in Section II. 

The growing availability of three-dimensional (3-D) volu- 
metric and spatio-temporal data increases the interest for 3-D 
two-channel designs. The dilation matrix associated to the 3-D 
face-centered cubic (FCC) sublattice is often proposed as the 
natural extension of the 2-D quincunx case [7]-[9]. However, 
we point out a fundamental difference between the 2-D quin- 
cunx and the 3-D FCC matrices: the 3-D FCC matrix does not 
correspond to a similarity transform (i.e., an angle-preserving 
transformation such as rotation or symmetry combined with di- 
lation) as does the 2-D quincunx matrix. Moreover, we show 
that, for more than two dimensions, it impossible to find a dila- 
tion matrix on the Cartesian lattice that is a similarity transform 
and at the same time leads to a two-channel filterbank. Con- 
sequently, when the 3-D FCC dilation matrix is cascaded, the 
Voronoi cells of the sublattices are not just obtained by a simple 
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rotation/symmetry and dilation. In particular, for the 3-D FCC 
case, we observe an alternation between three types of lattices: 
Cartesian, FCC, and quincunx in the x — z plane. This observa- 
tion can be useful to guide the orientation of the 3-D data before 
processing. These results are discussed in Section III. 

II. TWO-DIMENSIONAL QUINCUNX SUBSAMPLING 

We characterize the Cartesian lattice by its sites k = [fci 4-2] , 
ki, k2 & Z. An admissible 2x2 dilation matrix D, mapping k 
to Dk, must satisfy the following properties [10], [11]. 

1) The new lattice forms a sublattice of the Cartesian lattice; 
a trivial property when each element of D is integer. 

2) The magnitude of each eigenvalue A, of D must be 
strictly larger than 1 to ensure a dilation in each dimension. 

After applying the dilation matrix, the density of the laltices sites 
is reduced by a I'actor |<lc(, D|. The Voronoi cell of the lallicc is 
defined as the set of points closer to the origin than lo any other 
lattice site, and completely characterizes the lattice's geometric 
structure [12]. The dual or reciprocal lattice is characterized by 
the matrix that is the transpose of the inverse dilation matrix; the 
reciprocal lattice serves in the frequency domain to replicate the 
spectrum while its Voronoi cell can be considered as the Nyquist 
region [8], [13]. 

The quincunx sublattice consists of those lattices sites 
[ki k2]'^, where A;i -|- fc2 = 2n, n integer. There are two 
common choices for the corresponding dilation matrix: 



(1) 



Note that both matrices correspond to a similarity transform; 
i.e., a symmetry (Di) or a rotation (D2), combined with a 
dilation that reduces the number of sites by a factor of two 
(|detDi| = |detD2| = 2). For signal processing applications, 
an interesting property is to have the iterated sublattice coincide 
again with a dilated Cartesian lattice after a few iterations. This 
property also contributes to the construction of smooth wavelets 
[14]. For these examples, we have = 21 and D| = 161. 

The design of orthogonal filters for quincunx sampling is 
strongly determined by the shape of the Voronoi cell of the recip- 
rocal sampling lattice. More specifically, the goal is to design an 
orthogonal filter whose frequency response is reasonably close 
to the ideal filter characterized by the indicator fiinction of the 
Nyquist region of the sampling lattice in the frequency domain, 
which correspond to the Voronoi cell of the reciprocal lattice. 
Some filter and wavelet design examples can be found in [5], 
[7], [10], and [I5]-[19]. 
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(a) (b) 




Fig. 1 . Types of 3-D lattices encountered when the 3-D FCC dilation matrix is iterated. Left column: the Voronoi cell of the lattice in the spatial domain for the 
lattices (a) D^k, (c) k, (e) D^k (notice that the side in the x — z plane have length y/l,, while the sides in the y direction have length 2). Right column: the 
Voronoi cell of the reciprocal lattice in the frequency domain for the lattices (b) (D~'^)°k, (d) (D~^)^k, (f) (D~^)^k. (a) Cartesian lattice; (b) cartesian lattice; 
(c) FCC lattice; (d) BCC lattice; (e) quincunx lattice a; — z; (f) quincunx lattice — fx- 



Since the 2-D quincunx matrix corresponds to a similarity 
transform, its associated Voronoi cell in the spatial and the fre- 
quency domain is simply rotated or mirrored and dilated. There- 
fore, iterating Di and D2 alternates between a Cartesian and a 
quincunx sublattice. 

III. THREE-DIMENSIONAL FCC SUBSAMPLING 

The dilation matrix corresponding to the 3-D FCC sublattice 
is often considered to be the natural extension of the 2-D quin- 
cunx scheme. The corresponding most popular dilation matrix 
(cf. [4], [8], [10], and [20]) is 

[1 0 1] 
-1 -1 1 . (2) 
0 -1 oj 



It leads to a two-channel design (|detD| = 2) and coincides 
with a dilated Cartesian lattice after three iterations; i.e., we have 
= 21. 

The use of the FCC sampling scheme is motivated by the 
energy packing argument [21]. Indeed, together with the body- 
centered cubic (BCC) lattice, the FCC lattice has the densest 
lattice packing in 3-D: it retains a maximal energy proportion 
for signals with a spherical spectrum. 1 Therefore, as in the 2-D 
case, the shape of the frequency response of an orthogonally 
designed filter on the FCC lattice tends to the indicator function 

'In signal processing Literature, the term face-centered orthorhombic (FCO) is 
often used instead of FCC. However, the FCO lattice allows for a scahng along 
the coordinate axes and therefore the argument of maximal energy compaction 
for signals with a spherical spectrum is no longer valid. Since the scahng is 
rarely used, we prefer the use of the term FCC. 
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Properties of the 1 1 





spatial lattice 


reciprocal lattice sampling density 


3n 


Cartesian lattice 


Cartesian lattice 23«(1, 1, 1) 


3n + 1 


FCC lattice 


BCC lattice 23"( -?/2, ^) 


3n + 2 


Quincunx (x-z) 


Quincunx (fx-fz) 2^" ( v^, 2 , ^/2) 



with rn integer, wliicli is only possible for JV = 2. ■ 
Now we know that no 3-D dilation matrix with determinant 2 
exists that corresponds to a similarity transform. Nevertheless, 
it is interesting to take a closer look at the popular 3-D FCC 
matrix (2) and the consequences of not having this property. 
For this purpose, we determine the Voronoi cell at each itera- 
tion for both the spatial and the dual lattice [12], [21]. Since we 
have D'^ = 21, the iterated scheme alternates between three lat- 
tices, but due to the lack of being a similarity transform, with 
fundamentally different Voronoi cells. For D*^", n e N, we 
have the Cartesian lattice and the corresponding cubic Voronoi 
cell in the spatial and the frequency domain, see Fig. 1(a) and 
(b). For D"*""*"^, we have the FCC lattice in the spatial domain 
[Fig. 1(c)] where the Voronoi cell is a rhombic dodecahedron, 
and the BCC lattice in the frequency domain [Fig. 1(d)] where 
the Voronoi cell is a truncated octahedron. Finally, for D^""*"^, 
we obtain a quincunx sublattice in the a; — z plane, while the 
y dimension is aheady subsampled by 2; see Fig. 1(e) and (f), 
respectively. Logically, the Voronoi cell of the spatial domain 
gets larger when Ihc lallicc gets coarser, while the inverse hap- 
pens for the dual lallicc. The lalliccs arc also illuslrated by 2-D 
cuts in Fig. 2. The alternation between those lattices is summa- 
rized in Table I and has two important consequences. First, the 
optimal energy packing argument for orthogonal filter design is 
only valid for one of the three lattices. Second, the quincunx 
arrangement for each third lattice D'^""'"^ has a less isotropic 
sampling density, which also influences orthogonally designed 
filters. Applications can use this property to properly adapt the 
orientation of 3-D data before treatment; e.g., for spatio-tem- 
poral data, the temporal direction could be put along the y-axis. 



of the Voronoi cell of the BCC lattice, which is a space-filling 
polyhedron that fills the space in an optimal way for signals with 
a spherical spectrum. Some examples of 3-D filter design can be 
found in [10], [17], [19], [20], and [22]-[24]. 

However, there is a fundamental difference between the 2-D 
quincunx and the 3-D FCC dilation matrix. The latter does not 
correspond to a similarity transform. More generally, we can 
show the following. 

Proposition 1: For dimensions greater than 2, no dilation 
matrix can be a similarity matrix with determinant 2. 

Proof: Assume that D is an admissable dilation matrix in 
A'' dimensions that corresponds to a similarity transform. Then 
it should satisfy 

D^D = ml, (3) 
where m is an integer. Taking the determinant of (3) gives 

|detD| = V^. (4) 
On the other hand, by the two-channel design constraint, we 
need 

|detD| = 2. (5) 

Satisfying both constraints (4) and (5) would require 

TO= VI (6) 



rv. Conclusion 

The 3-D FCC dilation matrix is often considered to be the 
natural extension of the 2-D quincunx case. However, the di- 
lation matrices have a different nature. In this paper, we have 
shown that the popular 3-D FCC dilation matrix, as opposed to 
the 2-D quincunx matrix, does not correspond to a similarity 
transform. Moreover, for more than two dimensions, it is im- 
possible to build a two-channel dilation matrix on the Cartesian 
lattice which corresponds to a similarity transform. Still, the 3-D 
FCC dilation matrix can be useful to treat 3-D volumetric or spa- 
tiotemporal data. However, one should be aware that iterating 
the dilation matrix produces an alternation between lattices with 
a different degree of isotropy; i.e., Cartesian, FCC, and quincunx 
lattices. 
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